#ifndef hWW_H
#define hWW_H
#include "utility.h"
#include <stdio.h>
#include <iostream>
#include "mt2_bisect.h"

using namespace std;
using namespace mt2_bisect;

//****************************************************************************
extern "C" long double mH_finder_old(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	long double mH, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}
	//cout<<pt<<endl;


	mH = (EmM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 - (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) + (EpM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM)))
			/ (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 - (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) - (pLmM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2))
			* (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 - (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) - (pLpM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 - (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4
			* (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) + (EmM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) + (EpM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM
			- sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) + (pLmM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4
			* EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)))))) / (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) + (pLpM * sqrt(-((-1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM))) * (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4
			* pLmM * pLpM)))))) / (1 + (4 * EpM * pLmM + 4 * EmM * pLpM - sqrt(pow(-4 * EpM * pLmM - 4 * EmM * pLpM, 2) - 4 * (4 * EmM * EpM - pow(mW, 2)) * (pow(mW, 2) + 4 * pLmM * pLpM))) / (2. * (pow(mW, 2) + 4 * pLmM * pLpM)));

	cout << mH << endl;
	//cout << sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)) << endl;
	//cout << "**************" << endl;


	return mH;
}


//****************************************************************************
extern "C" long double mH_simple(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double M, alpha, mW;
	TLorentzVector pplus, pminus;

	mW = 80.4;

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		if (part[vectorstoadd[vect]] -> PID < 0) {
			pplus = vec_part[vectorstoadd[vect]];
		}

		if (part[vectorstoadd[vect]] -> PID > 0) {
			pminus = vec_part[vectorstoadd[vect]];
		}
	}

	alpha = pminus.Pz() / pplus.Pz();

	M = mW * (1 + alpha) / sqrt(alpha);

	return M;
}

//****************************************************************************
extern "C" long double mH_finder(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming a trial mH, and minize over the checker function over the trial mass.
	double mH, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;
	//mH = 160;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}

	float minchecker;
	float result;
	result = 0;
	minchecker = 8000;

	for (mH = 175; mH < 206; mH += 1) {

		float betaGuess;

		betaGuess = (EmM - EpM + EmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + EpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5)) * pow(pLmM - pLpM + pLmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + pLpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5), -1);
		float checker;
		checker = -4 * pow(mH, -2) * pow(EmM + EpM - (pLmM + pLpM) * betaGuess, 2) * pow(-1 + pow(betaGuess, 2), -1);

		if ((checker - 1) * (checker - 1) < (minchecker - 1) * (minchecker - 1)) {
			minchecker = checker;
			result = mH;

		}

	}

	return result;
}

//****************************************************************************
extern "C" long double beta_guess(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	double mH, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;
	mH = 190 * 1.0;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	float betaGuess;

	betaGuess = (EmM - EpM + EmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + EpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5)) * pow(pLmM - pLpM + pLmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + pLpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5), -1);

	return betaGuess;
}

//****************************************************************************
extern "C" long double mHguess_finder(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	double mH, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 80.4;
	mH = 190 * 1.0;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	float betaGuess;

	betaGuess = (EmM - EpM + EmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + EpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5)) * pow(pLmM - pLpM + pLmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + pLpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5), -1);
	float checker;
	checker = -4 * pow(mH, -2) * pow(EmM + EpM - (pLmM + pLpM) * betaGuess, 2) * pow(-1 + pow(betaGuess, 2), -1);
	//cout << checker << endl;
	//cout << sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)) << endl;
	//cout << "**************" << endl;


	return checker;
}

//****************************************************************************
extern "C" long double mHguess_finder2(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	double mH, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;
	mH = 170 * 1.0;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	float betaGuess;

	betaGuess = (EmM - EpM + EmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + EpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5)) * pow(pLmM - pLpM + pLmM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5) + pLpM * pow(pow(mH, -2) * (pow(mH, 2) - 4 * pow(mW, 2)), 0.5), -1);
	float checker;
	checker = pow(EmM - EpM - pLmM * betaGuess + pLpM * betaGuess, 2) * pow(-2 * EpM * (pLmM + pLpM) * betaGuess + 2 * EmM * (EpM - (pLmM + pLpM) * betaGuess) + pow(EmM, 2) + pow(EpM, 2) - pow(mW, 2) + 2 * pLmM * pLpM * pow(betaGuess, 2) + pow(mW, 2) * pow(betaGuess, 2) + pow(pLmM, 2) * pow(betaGuess, 2) + pow(pLpM, 2) * pow(betaGuess, 2), -1);
	//cout << checker << endl;
	//cout << sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)) << endl;
	//cout << "**************" << endl;


	return checker;
}

//****************************************************************************
extern "C" long double beta_comp(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	long double beta, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		// check if particle or antiparticle
		if (part[vectorstoadd[vect]]->PID == -13) {
			EpM = vec_part[vectorstoadd[vect]].E();
			pLpM = vec_part[vectorstoadd[vect]].Pz();
		}
		if (part[vectorstoadd[vect]]->PID == 11) {
			EmM = vec_part[vectorstoadd[vect]].E();
			pLmM = vec_part[vectorstoadd[vect]].Pz();
		}

	}

	beta = (0.5 * (4. * EpM * pLmM + 4. * EmM * pLpM - 1. * sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)))) / (pow(mW, 2) + 4. * pLmM * pLpM);
	cout << "**************" << endl;
	cout << beta << endl;
	cout << sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)) << endl;

	if ((beta > 1) || (beta < -1)) {

		beta = 10;
	}

	if (pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM) < 0) {
		beta = -10;
	}

	long double res;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	res = sum.Beta();
	float beta_sign;
	beta_sign = sum.Pz() / sum.E();

	cout << beta << " " << res << endl;
	float result;
	result = beta - beta_sign;
	return result;
}

//****************************************************************************
extern "C" long double beta_finder(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the beta of the Higgs assuming an mH using the four-vectors of two leptons
	long double beta, EpM, pLpM, EmM, pLmM;
	float mW;
	mW = 81;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			// check if particle or antiparticle
			if (part[vectorstoadd[vect]]->PID < 0) {
				EpM = vec_part[vectorstoadd[vect]].E();
				pLpM = vec_part[vectorstoadd[vect]].Pz();
			}
			if (part[vectorstoadd[vect]]->PID > 0) {
				EmM = vec_part[vectorstoadd[vect]].E();
				pLmM = vec_part[vectorstoadd[vect]].Pz();
			}

		}

	} else {
		cout << " this beta finder is meant only for two particles input" << endl;
		exit(1);
	}
	//cout<<pt<<endl;

	beta = (0.5 * (4. * EpM * pLmM + 4. * EmM * pLpM - 1. * sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)))) / (pow(mW, 2) + 4. * pLmM * pLpM);
	cout << beta << endl;
	cout << sqrt(pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM)) << endl;
	cout << "**************" << endl;

	if ((beta > 1) || (beta < -1)) {

		beta = 2;
	}

	if (pow(-4. * EpM * pLmM - 4. * EmM * pLpM, 2) - 4. * (4. * EmM * EpM - 1. * pow(mW, 2)) * (pow(mW, 2) + 4. * pLmM * pLpM) < 0) {
		beta = -2;
	}

	return beta;
}

//****************************************************************************
extern "C" long double y_positron_cm(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double rapidityofn, Ep, pLp, Em, pLm, E, pZ, Yp, Ym, YpCM, YmCM;

	TLorentzVector sum;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		if (part[vectorstoadd[vect]]->PID == 11) {
			//puts("trovato positivo");

			Em = vec_part[vectorstoadd[vect]].E();
			pLm = vec_part[vectorstoadd[vect]].Pz();

			Ym = 0.5 * log((Em + pLm) / (Em - pLm));
		}
		if (part[vectorstoadd[vect]]->PID == -13) {
			//puts("trovato negativo");

			Ep = vec_part[vectorstoadd[vect]].E();
			pLp = vec_part[vectorstoadd[vect]].Pz();

			Yp = 0.5 * log((Ep + pLp) / (Ep - pLp));
		}

		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	E = sum.E();
	pZ = sum.Pz();

	rapidityofn = 0.5 * log((E + pZ) / (E - pZ));
	//cout << Ym << endl;

	YpCM = Yp + rapidityofn;
	YmCM = Ym + rapidityofn;

	return YpCM;

}

//****************************************************************************
extern "C" long double rapidity_charge_asym(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************

	long double rapidityp, rapiditym, Ep, Em, pLp, pLm, result;
	//int vect = 0;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		//cout<<part[vectorstoadd[vect]]->PID<<endl;
		if (part[vectorstoadd[vect]]->PID > 0) {
			//puts("trovato positivo");
			Em = vec_part[vectorstoadd[vect]].E();
			pLm = vec_part[vectorstoadd[vect]].Pz();
		}
		if (part[vectorstoadd[vect]]->PID < 0) {
			//puts("trovato negativo");
			Ep = vec_part[vectorstoadd[vect]].E();
			pLp = vec_part[vectorstoadd[vect]].Pz();
		}

	}
	//cout<<pLm<<" "<<Em<<endl;
	rapidityp = 0.5 * log((Ep + pLp) / (Ep - pLp));
	//cout<<rapidityp<<endl;
	rapiditym = 0.5 * log((Em + pLm) / (Em - pLm));
	result = (rapidityp - rapiditym) / (rapiditym + rapidityp);
	return result;

}

//****************************************************************************
extern "C" long double delta_y_PmM(vector<int> vectorstoadd, TLorentzVector vec_part[], TRootLHEFParticle* part[]) {
	//****************************************************************************

	long double rapidityp, rapiditym, Ep, Em, pLp, pLm, result;
	//int vect = 0;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		//cout<<part[vectorstoadd[vect]]->PID<<endl;
		if (part[vectorstoadd[vect]]->PID > 0) {
			//puts("trovato positivo");
			Em = vec_part[vectorstoadd[vect]].E();
			pLm = vec_part[vectorstoadd[vect]].Pz();
		}
		if (part[vectorstoadd[vect]]->PID < 0) {
			//puts("trovato negativo");
			Ep = vec_part[vectorstoadd[vect]].E();
			pLp = vec_part[vectorstoadd[vect]].Pz();
		}

	}
	//cout<<pLm<<" "<<Em<<endl;
	rapidityp = 0.5 * log((Ep + pLp) / (Ep - pLp));
	//cout<<rapidityp<<endl;
	rapiditym = 0.5 * log((Em + pLm) / (Em - pLm));
	result = fabs(rapidityp) - fabs(rapiditym);
	return result;

}

#endif
